function grad = Gradient(u, boundary_condition, differencing_scheme, h, sigma)
% grad = Gradient(u, boundary_condition, differencing_scheme, h, sigma)
% 
% Gradient calculates the gradient of a function u : R^2 -> R (a 2D matrix)
%
% input:    - u : 2D image with scaler values
%           - boundary_condition    (DEFAULT = 'constant')
%           - differencing_scheme   (DEFAULT = 'central')
%           - h                     (DEFAULT = 1)
%           - sigma                 (DEFAULT = 0)
%
% output:   - grad : a 3D matrix, the kth partial derivative at i, j is stored
%             in grad(i,j,k)
%
% See also: Partial_derivative_x
%
% 2013-07-25 --------------------------------------------------------------
% 
% Matthias J. Ehrhardt
% CMIC, University College London, UK
% 
% matthias.ehrhardt.11@ucl.ac.uk
% 
% -------------------------------------------------------------------------

    if nargin < 2, boundary_condition = 'constant'; warning('boundary_condition not chosen. DEFAULT = %s', boundary_condition); end;    
    if nargin < 3, differencing_scheme = 'central'; warning('differencing_scheme not chosen. DEFAULT = %s', differencing_scheme); end;
    if nargin < 4, h = 1; warning('h for discretization not chosen. DEFAULT = %3.2e', h); end;
    if nargin < 5, sigma = 0; warning('sigma for smoothing not chosen. DEFAULT = %3.2e', sigma); end;

    grad = zeros([size(u) 2]);
    grad(:,:,1) = Partial_derivative_x(u, boundary_condition, differencing_scheme, h, sigma);
    grad(:,:,2) = Partial_derivative_x(u', boundary_condition, differencing_scheme, h, sigma)';
end